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We study probability distributions ol waves of topplings in the Bak-Tang-Wiesenfeld model on hy- 
percubic lattices for dimensions D > 2. Waves represent relaxation processes which do not contain 
multiple toppling events. We investigate bulk and boundary waves by means of their correspondence 
to spanning trees, and by extensive numerical simulations. While the scaling behavior of avalanches 
is complex and usually not governed by simple scaling laws, we show that the probabihty distribu- 
tions for waves display clear power law asymptotic behavior in perfect agreement with the analytical 
predictions. Critical exponents are obtained for the distributions of radius, area, and duration, of 
bulk and boundary waves. Relations between them and fractal dimensions of waves are derived. 
We confirm that the upper critical dimension Du of the model is 4, and calculate logarithmic cor- 
rections to the scaling behavior of waves in D = 4. In addition we present analytical estimates for 
bulk avalanches in dimensions D > 4 and simulation data for avalanches in D < 3. For D = 2 they 
seem not easy to interpret. 
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I. INTRODUCTION 

The sandpile model was introduced by Bak, Tang, and 
Wiesenfeld (BTW) [Q as a simple example of a slowly 
driven dissipative system exhibiting self-organized criti- 
cality (SOC). Although today many systems with SOC 
are known, it is considered as the prototype of such mod- 
els, and there is a huge literature devoted to it. Its 
theoretical understanding is crucially related to Dhar's 
discovery of its Abelian structure ||^ which allows exact 
calculation of many of its properties However, a 

complete analytical determination of the scaling behav- 
ior of avalanches is still lacking. Several approximation 
schemes, including a random walk approach Q , diffusion- 
like analogy [|| , renormalization group and a graph 
theory method were proposed, but led to different re- 
sults. In addition, computer simulations - which first had 
suggested simple scaling behavior together with standard 
finite-size scaling (FSS) - provide increasing evidence 
that the avalanche statistics is much more complicated. 
While most recent authors agree upon a breakdown of 
FFS, the detailed interpretation of their data is highly 
controversial among different groups [pO|-p^ . 

The standard FSS ansatz implies an asymptotic form 



(1) 



for the distribution of the number a of toppled sites in an 
avalanche (in other words, its "area"), where L is the size 
of the lattice, p{z) is a universal function, and Ta and Va 
are critical exponents. This ansatz implies simple scal- 
ing of all moments of a, (a") L"" with Un ctq + rti^a- 
Similar ansatzes should, according to this view, hold for 
the number of topplings s (which differs from a because 
sites can topple more than once in an avalanche) and for 
the radius and duration of avalanches. But recent inves- 
tigations |lO|-p^ show that the two-dimensional BTW 
model may be characterized by a multifractal behavior 
where different moments of a are governed by exponents 
cr„ which are not linear in n and are indeed not related to 



each other for different n. Different reasons for this have 
been proposed in ||l^,|ll| and in Notice that multi- 
fractality of avalanches can be proven for certain variants 
of the 1-d sandpile model |14|. 

Deviations from pure power laws had been seen already 
in early simulations, but were usually interpreted as fi- 
nite size effects due to avalanches which touch the bound- 
ary of the lattice. To illustrate that this is most likely 
not true, and that there is a real problem with simple 
scaling, we show in Fig. |l| the ratio Vsix, L)/Va{x, L) of 
the integrated distributions Vs{x,L) = dx' Ps{x' , L) 
ioi D — 2 and for different values of L. In these simula- 
tions, cylindrical boundary conditions were used (open at 
y — zLL/2 and periodic at a; = 0, L), and data were col- 
lected only for avalanches starting at y = 0. In this way 
we hope to have minimized boundary effects. Also, since 
we do not make separate fits to Ps{x, L) and Pa{x, L), we 
have none of the uncertainties inherent in such fits. Due 
to Eq. (|^) we would expect this ratio to scale as x'^'^^^" for 
X ^ minlL*^" , L"^' } « . According to analytical pre- 
diction ||l^ and recent large scale simulations p^ , p^ , p^ , 
the difference Ta — Tg should be in the range 0.024 to 0.08. 
The behavior seen in Fig. ^is rather different. Although 
the curves for different large L perfectly superimpose in 
a wide range, they are in this range not straight at all (as 
expected for a power law) , and their average slope in this 
universal range is much smaller. Very small differences 
Ta — Ts have been seen in several simulations using small 
lattices 18 . But it still disagrees with our data show- 
ing a lack of scaling even for avalanches which do not 
reach the boundary of the lattice. Most other variables 
show similar deviations from pure power laws in D = 2 
when scrutinized closely. 

In principle, one can expect that these deviations of 
scaling can be explained by assuming that the avalanche 
boundary advances like a pinned surface in a random 
medium. Unfortunately, this interpretation seems unten- 
able. As shown in ||l^ (see also [^), avalanches proceed 
in distinct waves of topplings. In each wave, any site top- 
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FIG. 1. Ratio 'Ps{x,L)/'Pa{x,L) of the integrated s- and 
a-distributions for two-dimensional sandpiles. According to 
the generally accepted FSS ansatz this should be a power law 
with exponent Ta — Ts in the range 0.024 to 0.08 in the region 
where it is independent of L. The dashed line is x'^"^. 



pies at most once. In the original version of the model 
waves overlap in time, but they can be disentangled by 
a simple trick so that at any time only a single wave 
propagates. Therefore, if at all, the arguments associ- 
ated with pinning effects should not apply to boundaries 
of avalanches but to the propagation of wave boundaries, 
and they would suggest that waves show complex behav- 
ior (notice that boundaries of successive waves are not 
simply related to each other 

But waves do behave simply, and do show simple scal- 
ing behavior. This is indeed the main message of the 
present paper. Our results extend analytical results 
derived in |19 2^,|2^ and large scale simulations made 
in 1^,0 . While the behavior of avalanches is complex 
and badly understood when typical avalanches are com- 
posed of many waves (which is the case for D = 2 and, 
to a much less degree, for D = 3), the behavior of single 
waves is simple and well understood. 

In particular, we show in the present paper that FSS 
works for waves of topplings. Since boundary avalanches 
(i.e. avalanches which start from an unstable boundary 
site) always consist of single waves Q , it applies also to 
them. Using the spanning tree representation of waves, 
the equivalence between spanning trees and loop-erased 
random walks, and rigorous estimations for the latter, we 
determine the critical exponents of their probability dis- 
tributions for all dimensions D > 2. We also use the wave 
statistics to confirm recent numerical [p5| and analyti- 
cal [pp^] predictions for the upper critical dimension of 
the BTW sandpile model. The upper and lower bounds 
for logarithmic corrections to scaling for four-dimensional 
waves are determined analytically and confirmed numer- 
ically. 

We discuss the possibility for investigation of avalanche 
distributions of the BTW model using the results ob- 
tained for waves. One of the key characteristics in this 
study is the average number {n)a of waves in avalanches 
of a given area a. Our simulations show that in D = 2 
this number grows, most probably, not slower than a 
power law with the exponent 0.17 (the first conjecture 



of was 1/6, and an analytical prediction of Q is 
1/4). This means that multiple topphngs substantially 
change the scaling behavior of avalanches in comparison 
to waves, which was indicated in many previous studies 
of the two-dimensional BTW model. 

In Z? = 3 the fraction of avalanches containing more 
than one wave is much less than in two dimensions. All 
accurate numerical estimations of the exponent lead to 
Ta — 1.33 which coincides with the exact exponent 4/3 for 
the wave distribution, which could mean that the aver- 
aged number of waves in an avalanche in _D = 3 grows not 
faster than logarithmically. On the other hand, consid- 
ering the numerical estimation of this number we cannot 
exclude its slow polynomial growth with the avalanche 
size. Then, the scaling behavior of avalanches could be 
corrected for the multiple topplings in large events, sim- 
ilar to the case D — 2. 

Finally, for D > A the upper logarithmic bound for the 
averaged number of waves in an avalanche [ p6[ implies 
that the distributions of avalanches obey asymptotic be- 
havior with the same exponents as for waves. 

This paper is organized as follows. In section || we 
remind basic definitions of the BTW model. Section |l| 
is devoted to detailed explanation of the construction of 
waves and their spanning tree representation. In sec- 
tion IV we derive the critical exponents of wave distribu- 
tions. In section ^ we discuss analytical results for the 
dynamical exponent and fractal dimension of waves. Re- 
sults of computer simulations are presented in section [v| . 



II. THE BTW MODEL 

We consider the _D-dimensional BTW model on a hy- 
percubic lattice of linear size L in which integer variables 
Zi > represent local energies. One perturbs the system 
by adding particles at randomly chosen sites according 
to 



+ 1. 



(2) 



A site is called unstable if the corresponding energy Zi 
exceeds the critical value 2D. An unstable site relaxes, 
its energy is decreased by 2D and the energy of the 2D 
nearest neighbors (nn) is increased by one: 



2D 



1. 



(3) 



(4) 



In this way, the neighboring sites may be activated and 
an avalanche of relaxations may proceed. If a boundary 
site topples, one or more particles leave the system. The 
avalanche of relaxations stops when all sites are stable 
again. 

One can introduce in a natural way different kinds of 
sub-avalanches, e.g. clusters of sites toppled not less than 
a given number of times or waves of topplings |jl^ . A 
relaxation event (an avalanche or sub-avalanche) is char- 
acterized by its size s (total number of topplings), area 
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a (number of distinct toppled sites), duration t (num- 
ber of parallel update steps until stable configuration is 
reached), and its radius r (e.g. the radius of gyration or 
the maximal distance between the origin and a toppled 
site) . The basic hypothesis of Bak et al. claimed that 
in the self-organized critical state the probability distri- 
butions of values s,a,t,r exhibit power-law behavior for 
intermediate values of x, 



(5) 



with a; € {s, a, i, r}. 

As we have seen, this hypothesis might not be true 
for complete avalanches, but as we shall see it does hold 
for waves. Scaling relations for the exponents Ts,Ta,Tt, 
and Tr can be obtained if one assumes that size, area, 
duration and radius of "typical" events scale as powers 
of each other, for instance 



t 



(6) 



Then the transformation law of probability distributions 
Pt{t)dt — Pr{r)dr leads to the scaling relation 



1 



7tr 



(7) 



Again we should warn that there is a crucial assump- 
tion underlying these relations, namely that conditional 
distributions Px{x\y) are narrow, and therefore Eq. (|^) 
holds with small deviations for most events. It was pro- 
posed in |12| that this might not be justified in D = 2, 



and this is indeed the main source of problems of this ap- 
proach. Let us ignore this for the moment and proceed 
nevertheless. 

The scaling exponents jxx' are important for the de- 
scription of the extent of avalanches and their propaga- 
tion. For instance the exponent indicates if multi- 
ple toppling events are relevant (7^0 > 1) or irrelevant 
(7sa = !)• The exponent jar relating the avalanche area 
to its radius r equals the fractal dimension D/ of the 
avalanche. Finally, the exponent "ftr is usually identified 
with the dynamical exponent z. 

If Eqs. (|[]^) are applied to waves, one has of course 
7sa = 1, but Jar and "ftr are non-trivial. Our main result 
states that Eqs. (H]^) do indeed apply to waves, together 
with the FSS ansatz Eq. (|l]). 



III. WAVES OF TOPPLINGS AND THEIR 
SPANNING TREE REPRESENTATION 

Dhar proved 1^ that all stable configurations can be 
classified as either transient or recurrent. The former 
can occur only during an initial transient period, but are 
irrelevant for the infinite time dynamics. He also for- 
mulated the so-called " burning algorithm" which, on the 
one hand, allows one to distinguish the recurrent states 
among all stable configurations, and, on the other hand, 
can be used for constructing a spanning tree representa- 
tion of any given recurrent state. According to this al- 
gorithm, which also proceeds in discrete time steps, any 



site i is "burnt" at time t if its energy Zi is larger than the 
number of its unburnt nearest neighbors at time t—1. In 
a stable configuration, only some of boundary sites can 
satisfy this condition at the first step, they can be inter- 
preted as origins of "fire" . Then the "fire" propagates if 
new sites become burnable at the second step. In this 
way, we burn the sites step by step, until no more sites 
can be burnt. If all sites of the lattice are burnt, the ini- 
tial configuration of energies is recurrent, otherwise it it 
transient. 

In order to obtain the spanning tree representation of 
recurrent configurations |15[| , each burnt site i is con- 
nected by a bond to one of the sites which had "set it 
afire" , i.e. which had caused its burning by burning itself. 
If there is more than one such site, one uses an arbitrary 
but definite set of rules where to place this bond. In ad- 
dition, one introduces a new site □ ( "sink" ) and connects 
it to all boundary sites. On these connections, bonds are 
placed to those sites which burn at i = 1. Then we can 
imagine the entire process to start by burning the site 
□ at time t = 0, and generating a rooted tree with root 
at □. If the state is recurrent, this tree spans the entire 
lattice. 

Majumdar and Dhar [|5| also noticed that the condi- 
tion for "toppling" of a site is essentially the same as the 
condition for "burning" : At each step, the site i topples if 
its energy Zi is larger than the number of those of its near- 
est neighbors which had not toppled in the step before. 
Therefore, the burning of a recurrent state is equivalent 
to a topphng process initiated "from the boundary" . It 
implies that if we add one particle to every boundary site 
(two particles on each corner, etc.), each site will topple 
exactly once during the ensuing avalanche. 

The burning algorithm gives a one-to-one correspon- 
dence between recurrent states and spanning trees. This 
allows one to calculate the total number of recurrent 
configurations, the energy probabilities and the energy- 
energy correlation functions p|-p|j27[|. 

The spanning tree representation can be constructed 
also for a certain class of unstable configurations ap- 
pearing during an avalanche. It was shown in ||l9| that 
avalanches in the BTW model can be decomposed into 
so called 'waves of topplings'. According to this con- 
struction, an avalanche is considered as a superposition 
of successive sub- avalanches. After perturbing the sys- 
tem at a given lattice site i, one allows it to relax, but 
prevents the site i temporarily from toppling a second 
time. After this first 'wave' all sites are again stable ex- 
cept, possibly, the site i. If i is unstable, a second wave 
is initiated by toppling it again. But a possible third 
toppling is again delayed until this wave is finished, and 
when it finally occurs, it triggers the third wave. The 
procedure is repeated until the site i is stable. Note that 
if the site i is on the boundary, the avalanche stops after 
first relaxation and consists of only one wave. More gen- 
erally, if the distance of i to the boundary is d, then any 
avalanche starting at i can have at most d + I waves. 

To obtain the tree representation of a configuration fol- 
lowing a wave which had started at site i, we introduce 
an auxiliary BTW model on a new lattice. In this lattice 
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FIG. 2. Spanning 2-tree representation of a wave. Toppled 
sites are marked by heavy dots. The origin of the wave is 
marked by a circle. The dotted lines indicate the boundary 
of the system. The latter is considered as a single additional 
site □, so that all non-toppled sites form a single connected 
tree. 



we connect i to the sink □. Since the site i in the new 
model is on the boundary, each avalanche starting at this 
site will consist of a single wave. Each avalanche in the 
auxiliary model corresponds indeed to a wave of some 
avalanche in the original model. Applying the burning 
algorithm, we can construct a spanning tree on the aux- 
iliary lattice representing the recurrent state of the new 
model. During this burning process, some branches of 
the fire will be independent of site i, but one branch will 
first pass from □ to i and then propagate further. It is 
this latter branch which coincides with the last wave of 
topplings in the original BTW model. 

Removing the bond between i and □ we obtain two 
trees, one having the root at i and the second at the sink 
□ . The first tree represents the wave and the second one 
corresponds to the sites not toppled in this wave. The 
tree with the root □ contains information about the con- 
figuration of stable sites not affected by the wave. We will 
call this union of two trees which cover the entire lattice, 
a spanning two-component tree (or, simply, a spanning 
2-tree) . According to this tree representation, waves with 
different configurations of untoppled sites are counted as 
different. An example of a spanning 2-tree is shown in 
Fig. |. 

The rigorous proof of the above construction is given 
in A similar decomposition of avalanches into "in- 
verse avalanches" was proposed by Dhar and Manna . 

An important fact concerning the wave statistics 
should be noted. Since all recurrent states of the auxil- 
iary BTW model have equal probability of occurrence , 
all waves in the original BTW model are also equally 
likely. 



IV. CRITICAL EXPONENTS AND GREEN 
FUNCTIONS 



Using the graph representation of waves, we can ex- 
press their probability distribution by the lattice Green 
function. Consider avalanches initiated by adding a par- 
ticle at the site i and spanning 2-trees with roots at i and 
□ representing waves of these avalanches. It was proven 
in ||l9| that the Green function dj is related to the num- 
ber iV(n)(y) of spanning 2-trees where the site j is in the 
same component as i: 



^(□)(y) 



N, 



(□) 



(8) 



where A^(n) is the total number of spanning trees with 
the root at the site □. 

Consider avalanches initiated by adding a particle at 
the site i. Since every recurrent state together with the 
perturbed site i completely defines the relaxation process, 
the number of different possible avalanches N^°'^ started 
at fixed point i equals to the number of recurrent states 
(or the number of spanning trees iV(n)). The number 
N^^^s^ of different waves started at the site i and covering 
the site j corresponds to the number iV(n)(y) of spanning 
2-trees having sites i, j on one of its components. We can 
therefore rewrite Eq. (0) as 



w) 



N, 



(a) ■ 



(9) 




is another formulation of the known result of 
that the expected number of topplings at site 
j due to adding a particle at the site i is given by the 
Green function Gy . 

Due to uniformness of the wave statistics mentioned at 
the end of section III, the probability that a wave W{i) 
starting at the site i covers the site j is equal to the 
fraction of waves having this property: 



P{j e W{i)) = 



(10) 



where N^^™^'' denotes the total number of waves starting 



at i. 



Combining Eqs. @ and (|T^), we write 



P{j e W{i)) = 



G^ 
Gi 



G{r) 
G(0)- 



(11) 



where we use the notation G(r) for the Green function 
Gij if the points i and j are separated by the distance r. 
On the other hand, this probability can be represented 

as 



/>oo 

PijeWii))^ / P'-^^R) pB.{r) dR 

J r 



(12) 



where P*"''(i?) is the probability that the linear extent 
of a wave is R, and PR{r) denotes the density of sites 
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covered by such a wave. The density pR{r) tends to 1 
for large R if waves are compact and isotropic, and is 
a function of r if waves are fractal. Asymptotically for 
r, it varies as 



^(boundary) 



D 



(13) 



where d/ is the fractal dimension of waves and D is the 
Euclidean dimension of the lattice. 

Suppose that the probability distribution of the wave 
radius r has a power law asymptotics similar to that for 
avalanches [Eq. (H)] , 



(14) 



Then, the probability distribution P^"'>(V > r) scales 
with the exponent r^"> — 1. Using Eq. (|l3|) we get the 
asymptotic behavior of the probability in Ihs of Eq. (112) 



P{j e W{i)) - r- 



> + l+df-D 



(15) 



The asymptotics of the bulk Green function (see for 
instance [p8[) are given by 



G(r) 



Inr 



(16) 



for Z? 2 
for Z? > 2 

which reveals that the radius exponent t^™' for waves is 



=df~l. (17) 
Using Eq. (|^) we can derive the exponents of the wave 



area 



and duration 



= 2- 



= 1 



(18) 



(19) 



respectively. 

For avalanches started at a distance b from the bound- 
ary, we need the boundary Green functions which can be 
calculated by the method of images: 



ln|r + b| -ln|r-b| 



for D = 2 



G(r) 



|r-bp--D - |r + bp--D for D > 2 



(20) 



where b is the vector perpendicular to the boundary. On 
any "equipotential" surface G(r) = const characterized 
by a length scale ^ and a volume a ~ , this boundary 
Green function scales as 



l-D 



(21) 



If we now replace Eq. (|T^) (which is appropriate only 



-D 



for isotropic cases) by its generalization p ^ 
we arrive at the exponents for waves starting near the 
boundary: 



(boundary) 



(22) 



(23) 



(24) 



(here and in the following we use superscript (w) for bulk 
waves and (boundary) for boundary waves, and use sym- 
bols without superscripts for avalanches). We see that 
both the bulk and boundary wave exponents are deter- 
mined by the scaling exponents df and z. The dynam- 
ical exponent z can be related to the fractal dimension 
of the "chemical path" on a spanning tree jl^] which, in 
turn, is equivalent to the dimension of the loop-erased 
random walk (LERW) |^^. As to the fractal dimension 
of waves d/, it was proven for all dimensions that a set 
of untoppled sites, which are completely surrounded by 
toppled sites, corresponds to a forbidden subconfigura- 
tion ||l^,^. However, this fact does not prevent waves 
from being fractal. They still could display either non- 
fractal or fractal behavior depending on the dimension- 
ality D. In the next section, we will show that d/ is also 
closely related to properties of LERW, more precisely to 
the intersection probability between a LERW and a sim- 
ple random walk. 



V. LOOP ERASED RANDOM WALKS, 
DYNAMICAL EXPONENT AND FRACTAL 
DIMENSION OF WAVES 

In this section, we derive analytical estimates for crit- 
ical exponents of waves using their spanning tree repre- 
sentation and equivalence between a chemical path on a 
spanning tree and LERW. 

Consider an unrestricted random walk on a hypercu- 
bic lattice. The LERW introduced by Lawler is ob- 
tained from the simple random walk by deleting all loops 
along the path. The chemical path between two sites of 
a tree is the unique path along the tree edges connecting 
these sites. Majumdar |2^ has shown that the chemical 
path on a spanning tree is statistically equivalent to the 
LERW, i.e., the average distance r between the starting 
point and the position after I steps scales as r ~ with 
the same exponent for both of them. 

In Z? = 2, the exponent — 4/5 is known exactly 
from conformal field theory |31]. In _D = 3, mmicrical 
estimates yield i' w 0.616 pz 



33]. In D 



4 which is the 

upper critical dimension for the LERW, = 1/2 and the 
simple scaling law has logarithmic corrections For 
D > A, the scaling behavior of the LERW and chemical 
paths is given by the trivial value i' — 1/2, as the effects 
of self-intersections become negligible above the upper 
critical dimension. 

Returning to the BTW model, we notice that sites 
which topple at a given step of wave propagation coin- 
cide with sites deleted at the same step of the burning 
process, if it is started at the origin of the wave. Since 
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while a lower bound was obtained in [P6|: 




FIG. 3. A sketch of a wave initiated at the site i and con- 
taining the site j. The chemical path between the two sites 
on the tree representing the wave is r{i,j). The second path 
r(fc, □) corresponds to a random walk which starts at the site 
k and escapes the path T{i,j) until it is trapped at the sink. 



(Inr 



,1/2 



(lni?)i/3- 



(28) 



We can see that PR{r) approaches 1 when R ^ oo, r 
fixed. The only fractal dimension which is consistent 
with both upper and lower bounds Eqs. ( pTjp^ ) is 4, but 
there are logarithmic corrections. 

For D < A the lower bound Eq. (|28| ) becomes stronger 
because of increasing intersection probability Pint ikD\ij). 
Thus, we conclude that the fractal dimension of waves is 



D for D< 4 



for £> > 4 



(29) 



which means that the upper critical dimension for waves 
is 4. 

Therefore we can calculate from Eqs. HH) the 

exact values of all exponents for bulk and boundary waves 
for all dimensions D > 2, with a single exception. This 
exception is the exponent of the duration distribution in 
D — 3, for which we need the value of fLERw(^ = 3) 
which is not known exactly. 



the burning process generates a tree, there exist a unique 
path from the root i of the tree (the site where the wave 
is initiated) to one of the last toppled sites if oi the wave. 
The number of update steps is given by the number of 
edges in this path. Thus, the duration t of the wave 
equals to I, the length of the chemical path from i to if 
on the tree and, therefore, the dynamical exponent jtr of 
waves is given by z = 

In order to find the fractal dimension of waves df^ we 
use the proposition proved in |^^. For this we take a 
site k at distance r — \k — i\ from i and a site j at 
distance i? = |j — i| > r, together with some paths r(z,j) 
connecting i with j and r(A;, □) connecting k with the 
sink □ (see Fig. |^) . Then the density of sites at distance 
r from i, in waves of radius R > r starting at site i, is 
given by [|6|: 



PRir) - Pint{kn\ij) 



(25) 



where Piat{kD\ij) is the probability that r(fc,n) inter- 
sects the path T{i,j), averaged over all j, all paths from 
i to j, all k, and all T{k, □). 

Using the known estimations of the intersection prob- 
abilities we can obtain from Eq. (|5|) the following 
upper bounds. For Z? > 4, we have 



p{r) = lim pR{r) < Cir 

R^OG 



4-D 



(26) 



From Eq. (|13|), we can see that the fractal dimension of 
waves df < 4 for all D > 4. 

For D = 4, the upper bound reads 



ln(l + a) i?2 
Pfl(r)<C2 , a^ — 



(27) 



VI. COMPARISON WITH NUMERICAL 
SIMULATIONS 

A. D=2 

In this subsection we present the results of numerical 
simulations of bulk and boundary waves in 13 = 2. For 
these the standard FSS works well. For any of the ob- 
servables x — a,t, and r it can be written as 



(30) 



with (3x =Txi^x IQ- 

The functions gxiz) should be universal (i.e., does not 
depend on the type of lattice) . But for large values of z 
they do depend on the type of boundary conditions (open 
on all four sides or cylindrical, i.e. open in one direction 
and periodic in the other) and on the aspect ratio (square 
or rectangle with sides Li ^ L2). We verified that the 
exponents were independent of boundary conditions and 
aspect ratios. We verified also that all results were un- 
changed if we threw in the sand grains with non-uniform 
density, provided this density was everywhere non-zero. 
The latter was very useful since it allowed us to obtain 
much improved statistics from either the boundary or the 
central region. 

Since d/ = 2 for 13 = 2, one has Va = df = 2 and 
Vr = 1, and therefore = ^xr |3^- The results of the 
previous section, together with vt — z = I/i^lerw = 5/4, 
give 



2ri' 



ZTl 



5/4, 



(31) 
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FIG. 4. Finite-size scaling plot of the -wave area distribu- 
tion Pa"^[a) for bulk waves in D = 2. The perfect data col- 
lapse shows that d/ = 2 and /3^™' = 2, as predicted analyti- 
cally. The dashed line represents a power law with exponent 
^(vi) _ rpj^g factor InL comes from the normalization of the 
distributions. The inset verifies the scaling ansatz Eq. (111). 



for bulk "waves and 

^(boundary) _ 



2t( 



^(boundary) ^ ^^^(boundary ) ^ Q j 



(32) 



for boundary waves. 

The finite-size scaling plot for the area distribution 
P^™' (a) of bulk "waves in = 2 is sho"wn on Fig. ^. In the 
inset of this Figure, as "well as in the insets of plots for 
other distributions of waves, we show the collapses ac- 
cording to ansatz of Eq. (||). Taking /J^™' = 2 and df = 2 
we see a perfect data collapse. The finite-size scaling 
ansatz of the duration distribution is plotted in Fig. H. 



FIG. 6. Same as Fig. y, but for boundary waves in _D = 2. 
Again the data collapse is obtained with the predicted values 
^^boundary) ^ 3 d/ = 2. Thc dashcd line corresponds to 
^^boundary) _ 3^2^ ^s prcdlctcd thcoretically. 



These data confirm that waves are not fractal, and that 
their duration is as predicted by the correspondence with 
spanning trees and loop-erased random walks. 

As was mentioned above, the avalanches started at the 
boundary consist of a single wave, so for this type of 
avalanches the distribution of avalanches coincides with 
that of waves. The finite-size scaling plots for the area 
and duration distributions for boundary waves are shown 
in Fig. ^ and Fig. ^ respectively. Again all predictions 
are verified, in particular we see that df ~ 2, i.e. also 
boundary waves are not fractal in D = 2. 

Finally, let us consider avalanches starting at a fi- 
nite distance from the boundary. The crossover from 
the boundary to bulk behavior of the wave distribution 
should be described by Eq. (|20|). More precisely, the 
equipotential surfaces in D = 2 are circles Pq| with ra- 



10° 



10" 



-J 10 - Z 



10 



10° 




, -5/4 

L t 



FIG. 5. Finite-size scaling plot of the wave duration dis- 
tribution P/™' it) for bulk waves in D = 2. Here, the data 
collapse is achieved with /S^"^' — u^™^ — 5/4 which confirms 
again Eq. (pi|). The dashed line represents a power law with 
exponent = 1. Again, the normalization factor InL is 
needed for a good data collapse. 
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FIG. 8. Area distributions of waves initiated at different 
distances b from the boundary for a system of size 256 x 1024 
with cyhndrical boundary conditions. Each curve is averaged 
over approximately 10^ avalanches. The dashed lines corre- 
spond to the bulk (ri™' = 1) and boundary (T-^i»ou„dary) ^ g^g) 
scaling behavior, respectively. The distributions are multi- 
phed by h-^N(h),N{h) = j P^'"\a\h)da in order to have the 
curves collapsing for large a. Inset: the rescaled distributions 
for fe = 8, 16, 32, 64 demonstrate that the crossover from bulk 
to boundary behavior takes place at values of area of order 



dius ^ and G - ln[(6V^2 ^_ ^)i/2 _ 5/^]^ ^j^g scaling 
region where p(^) = 1, we have therefore P^"''(a|6) = 
— {da/d^)~^dG/d^ which gives 



hjo?!'^ for a > nb"^ 
1/a for a < nb'^ 



(33) 



To check this, we simulated the BTW model with cylin- 
drical boundary conditions, and collected data for waves 
started at distance b from the open boundary. The re- 
sults are plotted in Fig. ||. For small a we see indeed the 
bulk behavior which crosses over to the boundary behav- 
ior a~^^'^ for b^ < a < L^. In the latter region we also 
see clearly the linear dependence on b. 

The above shows that our understanding of waves in 
the two-dimensional BTW model is basically complete. 
In contrast, and in spite of numerous efforts, the scal- 
ing behavior of avalanches in the two-dimensional BTW 
model is still an open problem. This is due to multiple 
topplings. The average number of waves in an avalanche 
scales as [0 



Ini 



(34) 



There are also several results known about correlations 
in the sizes of successive waves |,|3,|^,|9|. Nevertheless, 
even most basic questions such as the distribution of n 
or the dependence of n on the area a are not yet solved. 

Equation (|4|) would be most easily explained if P„ {n) 
were simply Xjr?. Present data seem to agree with 
this for the largest lattices (Fig. H), but actually the data 
are better fited with a power Xjr?-^ than with Xjri} (see 
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FIG. 9. Probability distribution of the number of waves n 
in an avalanche for D — 2. The dashed line corresponds to 
the power- law behavior Pn{n) ~ as predicted in |l5|. The 
data were collected from avalanches initiated at the center 
region of the square lattice with cylindrical boundaries. The 
inset shows the same data multiplied by There, the dashed 
line is oc n~^'^ . 



inset). Similar results are obtained for the average 
number of waves in avalanches with fixed a (Fig. |lO[ ). 
Although they seem to scale like a power of a, as assumed 
in [l5|] , a closer study shows significant deviations which 
seem hard to explain as finite size effects. 

There are several recent papers which try 

to explain these problems by unexpected features of 
avalanches which reach the boundary. But data such 
as those shown in Figs. |l|,|l^ indicate that there are al- 
ready problems with avalanches which do not reach the 
boundary. 
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FIG. 10. Average number of waves (n)a as a function of 
the avalanche area a for D = 3 and various system sizes L. 
Although the main figure looks straight at first, the inset dis- 
playing the rescaled average shows significant deviations from 
the assumed pure power-law behavior |l^. The data were 
collected from non-dissipative avalanches initiated at the cen- 
ter region of the square lattice with open boundaries. Thus 
the curvature seen in the inset cannot come from avalanches 
reaching the boundary. 
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B. D=3 



Simulations of waves of the BTW model in 13 = 3 
also give good agreement with our analytical results 
of Sections IV and 0. 

and ri™' = 1 + i^ler 



For bulk waves 



we now have 
+ 1 « 2.623, 



1 + J^LERw ~ 1.616. For boundary waves, the 
corresponding numbers are ^^^oundary) _ ^^boundary) _ 
5/3, /^(•'-"'■^"y) = 1/z/i.ERw + 2 « 3.623, and ri"''""'"'^ = 
l + 2fLERw ~ 2.232. For these /3 values, the data collapses 
of bulk (Fig. |ll|,|l^) and boundary (Fig. [T^Jl^ ) waves are 
perfect. They confirm also the analytical predictions for 
the T exponents, verifying in particular that the waves 
have fractal dimension df — 3. 

Due to the rarity of multiple topplings, avalanche dis- 
tributions coincide within the displayed accuracy with 
wave distributions, when plotted as in Fig. ID and Fig.|l|. 
In order to show significant results for multiple topplings, 
we have to present the data differently. In Fig. |l5| we plot- 
ted the average number of waves at fixed a, {nja, against 
log a. Neither using a logarithmic scale for {n)a (main 
figure) nor a linear scale (inset) give perfectly straight 
lines. Thus the data can be interpreted either as a power 
law with a very small exponent. 



0.06 



(35) 



or as a logarithmic growth. 

In the latter case, we would of course have r^™' = Tq = 
Tg. In the opposite case of a power law with exponent 
a « 0.06 we can give crude estimates for the differences 
between these r exponents, using some heuristic assump- 
tions. 

The first assumption is that different waves in the same 
avalanche involve essentially the same sites. If this is 
true, we should have P^"'(a)da « {n)aPa{o-)da. Using 
this together with Eq. (35), we find Tq = 4/3 + a « 
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FIG. 12. Finite-size scaling plots of the duration distribu- 
tion Pj™\t) for bulk avalanches in _D = 3. Using compact 
avalanches {z = 1/i^lerw) we get good data collapses and the 



resulting Tj exponent agrees with the theoretical prediction 
(dashed lines). 



1.39. Since the basic assumption here is most likely not 
justified, this is only a very crude (and most likely too 
large, in particular since the growth of {n)aPa{a)da could 
be logarithmic) estimate for the difference between Tq 
and r^™' . 

An estimate for the difference between and is 
obtained as follows. An upper bound for the size of 
an avalanche of area a containing n waves is s < na. 
Therefore, the assumption s ~ {n)aCL leads to the maxi- 
mal difference between the area and size exponents. Us- 
ing Eq. (H) we get a > sV(i+"). Then, Eq. (g) gives 
(ts — 1) > (tq — 1)/(1 + a) and we can conclude that 



-(Ta- 1) < 0.02. 



(36) 
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sulting Ta™' exponent agrees with the theoretical prediction 
(dashed lines). 
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C. D=4 

As the dimension of the BTW model increases, multi- 
ple toppling events in avalanches become more and more 
rare. For D > 4 it was shown in that {n)a grows 
not faster than logarithmically, i.e. a = 0. As we al- 
ready mentioned in the previous subsection, this means 
r^™' = Ta = Ts, i.e., the scaling behavior of waves and 
avalanches in D = 4 is characterized by the same expo- 
nents and scaling functions. 

At the upper critical dimension = 4 logarithmic 
corrections to the simple scaling are essential. The prob- 
ability distributions of the radius, duration, area and the 
scaling relations between them are expected to have the 
form (cf. p5|) 



Prir) 



(In r) 



(Ini)''* 



Paia) 



(Ina)^" 



j3/2 



(37) 



and 



A direct verification of such slight differences between 
the T exponents could be tried by plotting ratios of the 
distributions, as was done in Fig. |l| for I? = 2. We do 
not show any such ratio here, since they are all very close 
to 1 in the region where the distributions should follow 
power laws, and the deviations from 1 seem not to be 
simple powers. 
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FIG. 15. Average number of "waves {n)a as a function of 
the avalanche area a for D — 3 and various values of the sys- 
tem size L. In order to minimize finite size effects, cylindrical 
boundary conditions "with one open and t"wo periodic direc- 
tions "were used and one half of all sand grains "were thro"wn 
onto the central planes y = L/2 and y — L/2 + 1. The 
y-axis is logarithmic in the main plot, and linear in the inset. 
Neither way of plotting gives perfectly straight lines in the 
region where that data for different L collapse. Although the 
main figure looks more straight at first sight, a more careful 
inspection shows a slight downward curvature. 



(Inr) 



Na 



t 



(Inr) 



Nt 



(38) 



respectively. 

The exponents of logarithmic corrections Xr, Xa, Xt, 
Na, and Nt obey the scaling relations Eq] 



(39) 



Xr = Xa+ Na/2, Xr = Xt + Nt 



which follow straightforwardly from Eqs. ( |37| ) and (|3J 

Using arguments similar to those at the end of the pre- 
vious subsection, we obtain an inequality for the expo- 
nents of logarithmic corrections for waves and avalanches 



Xn, ^ Xc 



(40) 



This allows us to compare below analytical estimations 
for waves with numerical results for avalanches. 

It follows from Lawler's results discussed in Sect. ^ 
that Nt — 1/3 exactly. But as with all logarithmic cor- 
rections, a numerical verification is not easy. The main 
reason is that the logarithms are never very much larger 
than 1, even for the largest simulations. Therefore the 
next-to-leading terms (which are typically suppressed by 
powers of the same logarithms) are in general not neg- 
ligible. In view of this, the disagreement with recent 
simulations which had suggested Nt « 1/2 should 
not be taken serious. Data for the mean squared radius 
of avalanches with fixed duration t are shown in Fig. 
More precisely, since we expect 



A^A^(lni)i/3, 



(41) 



we plotted [{r)^/t]'^ against Int. Apart from very large t 
when the finiteness of the lattice makes itself seen, we ob- 
serve essentially a straight line (which is a bit fortuitious 
since there are also 1/t corrections which are important 
for small t). At the same time, a power law dependence 
(r)t t^" with 1/ > 1/2, as would be expected if Dc > 4, 
seems ruled out. 



submitted to Phys. Rev. E 



11 



0.10 



0.08 - 



0.06 - 



0.04 - 



0.02 - 



0.00 



L=32,64,96,128,144 / 

/ — — ' 


D=4 


//V' \ 

— \ 









10 



100 

t 



1000 10000 




FIG. 16. The gyration radius of avalanches as a function 
of their duration. We plot the sixth power of the average 
rescaled gyration radius, [r^t"^]^, in a logarithmic diagram, 
since this should result in a straight line according to Eq. 
Such a linear regime is indeed observed (dashed line; its slope 
and intercept are not predicted by theory), and it increases 
with the system size L. 



For the other exponents Xr,Xa,Xt, Na we can only give 
inequahties from the analytical results of section The 
upper bound Eq. (|27| ) for the density of waves leads to 
the relation 



p{r) ^ (Inr) 



-<5 



S>1. 



Using this asymptotics and Eq. ( JTl] ) we get 



P^'M - G{r)/p{r) 



(42) 



(43) 



which gives Xr = 6- The area a of a wave scales in leading 
order as 

„4 



(44) 



giving Na = 5. Hence, from Eq. (^9|) we have 

Xa = 5/2 > 1/2. (45) 

In order to verify these predictions - and to verify, 
in the first place, that deviations from power laws with 
the mean field exponents = 3,Tf = 2, = 3/2 can- 
not be eliminated by changing these exponents - we per- 
formed extensive simulations. Numerical data of the size 
distribution Ps{s) are shown in Fig. 17, Lattice sizes 
range from L = 32 to L = 144. After multiplying with 
s^/^, we see indeed no indication that the remaining s- 
dependence follows a clean power law. In order to find 
the expected logarithmic corrections, we multiplied these 
data by (In s)^" , with several trial values for the exponent 
Xa- Taken at face value, this would give best fits with 
Xa ~ 0.25. In view of inequality ( p9| ) and of the diffi- 
culties in obtaining correct logarithmic corrections men- 
tioned above, we propose that indeed Xa = 1/2. From 
the relations Eq. ( p9| ) we get then Xr — Na = 1 and 
xt = 2/3. 



FIG. 17. Area distributions for avalanches in D = 4. In or- 
der to render the plots more significant, first of all the dom- 
inant s-dependence was removed by multiplying with s'^^^. 
Then, in order to check whether the remaining s-dependence 
in the scaling region is compatible with logarithmic correc- 
tions as proposed in Eq. (|3^) , we also divided by powers (In s)^ 
and shifted the resulting curves horizontally and vertically in 
order to avoid overlaps. The best agreement is found with 
X ^ 0.25. 



VII. CONCLUSIONS 



We have studied probability distributions of waves of 
topplings in the BTW model on D-dimensional rectan- 
gular lattices for D > 2. We have proved analytically 
that waves as well as boundary avalanches do exhibit 
critical behavior and that their probability distributions 
display power law asymptotics. We have derived exact 
values of critical exponents of these distributions. We 
have proven analytically that the upper critical dimen- 
sion of the BTW model is = 4, showing that pre- 
viously observed deviations from mean field behavior at 
Z? = 4 |17 3^ are due to logarithmic corrections. All 
these results have been confirmed by extensive numeri- 
cal simulations. During these simulations we have also 
verified that wave distributions follow the standard finite 
size scaling ansatz. The exponent of the leading logarith- 
mic correction to the distribution of avalanche life times 
(or, more precisely, lifetimes of waves) has been derived 
exactly from known asymptotics of loop-erased random 
walks. Estimations are given for the exponents of the 
logarithmic corrections to the other distributions. 



We therefore have now a rather complete picture of 
the dynamics of single waves in the BTW model for all 
dimensions. For D > 4 this means that we also under- 
stand avalanche dynamics, since multiple topplings are 
so rare that they can be neglected. For D = 2 the latter 
is certainly not true, and our understanding of avalanche 
dynamics is still incomplete. For D = 3, finally, multiple 
topplings represent a small but not negligible effect, and 
we have hope that the problem will be solved soon. 
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